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NEW ONE-STEP INTEGRATION METHODS OF HIGH-ORDER ACCURACY 

APPLIED TO 

SOME PROBLEMS IN CELESTIAL MECHANICS 


By 

Erwin Fehlberg 

SECTION I. INTRODUCTION 

1. The methods for the numerical integration of initial value problems in ordinary 
differential equations can be divided into two classes — multistep methods and 
one -step methods. 

2. Multistep methods were developed as early as the nineteenth century, mainly 
for astronomical problems. As their name indicates, these methods use the 
information from several backward (or also forward) computation steps in cal- 
culating the solution for the current step. Multistep methods (such as the 
ADAMS, GAUSS, COWELL, etc. methods) are very useful in problems for 
which the numerical integration can be performed with a constant step size. 
Since many such problems are encountered in astronomy, it is quite natural 
that astronomers have developed a number of powerful multistep methods . 
Moreover, since these multistep methods can be extended to any order of accu- 
racy (simply by adding higher-order difference terms to the formulas) and 
since, in general, they require only one or two evaluations of the differential 
equations per step, they seem to represent a rather rapid and economical 
integration procedure. 


3. 


However, multistep methods do have certain inconveniences and disadvantages: 



A. They are not self-starting but require a special starting procedure. 

B. Halving the step size during the computation requires time- 
consuming iterations to build a new difference scheme. 

C . The truncation error for multistep methods is larger than for one- 
step methods of corresponding order. Multistep methods therefore 
require a smaller step size than corresponding one-step methods. 

D. The classical multistep methods are, for stability reasons, of only 
a mediocre order of accuracy, considering the number of steps in- 
volved. Although a k-step formula for the solution of a first-order 
differential equation contains 2k + 1 constants, there exists no numer- 
ically stable k-step formula of an order exceeding k + 1 (for odd k) 

or k + 2 (for even k). This means that the stability requirement 
reduces the possible order 2k of such a k-step formula by k - 1 or 
k-2, respectively. 

4. Only recently, W. B. GRAGG and H. J. STETTER [14] have succeeded in elim- 
inating the stability restrictions of the classical multistep methods by introduc- 
ing into the formulas one extra non- step point. Such modified multi step 
formulas have been published by J. C. BUTCHER [ 8]. His paper contains 
numerically stable k-step predictor-corrector formulas of order 2k+ 1 for 

k ^ 6. Since our paper deals with one-step methods, we shall consider one of 
J. C. BUTCHER'S new multistep formulas in Appendix C for comparison only. 

5. In one-step methods, no information obtained from previously computed steps 
is required. Most one-step methods are of the RUNGE-KUTTA type. In 
RUNGE-KUTTA formulas the necessary information is obtained by repeated 
evaluation of the differential equation at intermediate points somewhere between 
the initial and the end point of the current step. 
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Since the standard RUNGE-KUTTA formulas are of fourth-order accuracy only— 
the truncation errors being of the fifth order of the step size — several exten- 
sions of these formulas to higher-order accuracy have been achieved in the 
last decade. We mention in this respect the work of A. HUTA [15], [16], 

J. C. BUTCHER [1], [2], [3], [6], and E. B. SHANKS [18], [19]. The last 
author has derived the most accurate RUNGE-KUTTA formula to date — an 
eighth-order formula based on 12 evaluations per step. 

6. Like all one-step methods, the RUNGE-KUTTA method is self-starting and 
the integration step size can be changed at any time and can immediately be 
accommodated to the local conditiohs of the problem under consideration. In 
this respect, RUNGE-KUTTA methods are well-suited to problems that re- 
quire frequent changes in the step size. However, RUNGE-KUTTA methods 
also have certain disadvantages. They are time-consuming, since they re- 
quire a relatively large number of evaluations per step of the differential equa- 
tions. Moreover, no economical method of step-size control seems to exist 
for RUNGE-KUTTA formulas. Apart from somewhat doubtful rule-of-thumb 
control procedures, there exists only L. F. RICHARDSON'S well-known method 
of the deferred approach to the limit. This method is quite reliable, but it 
doubles the computational effort merely for the benefit of step- size control. 

7. This paper will describe some one-step methods that the author has developed. 
They are essentially a combination of power series expansions and RUNGE- 
KUTTA methods. When applicable, our formulas have definite advantages 
compared with standard RUNGE-KUTTA formulas: they yield any order of 
accuracy one might desire; they require only very few evaluations of the 
differential equations; they include a very simple and economical method of 
step-size control. We shall describe these new RUNGE-KUTTA methods in 
detail in Section III. 
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8. In Section IV we shall derive the equations of motion of the restricted prob- 
lem of four bodies. Since this problem is of some practical interest in 
astronautics, in Sections V and VI we shall apply our new RUNGE-KUTTA 
procedure to this problem as well as to the restricted problem of three bodies . 

9. For comparison we shall report, in the appendices, our experience with some 
new methods developed by other authors. These appendices will give brief 
descriptions and applications of the following methods: 

A. The explicit RUNGE-KUTTA formulas of E. B. SHANKS (Appendix A). 

B. The implicit RUNGE-KUTTA formulas of J. C. BUTCHER (Appendix B). 

C. The modified multistep method of J. C. BUTCHER (Appendix C). 

10. A short abstract of a part of this paper was presented at the Congress of the 
International Federation of Information Processing (IFIP) in New York in 
May 1965 [13], This joint presentation, by S. FILIPPI and the author, also 
included work by the authors on the LIE series method. Dr. FILIPPI intends, 
at a later date, to publish the results of his LIE series investigations as a 
NASA Technical Note. 
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SECTION II. POWER SERIES EXPANSION METHOD 


11. The integration procedures that we shall describe in this and the following 
section are based on a power series expansion for the solution of the differ- 
ential equations under consideration. Such a power series expansion requires 
a repeated total differentiation of the differential equations with respect to 
the independent variable in order to obtain the coefficients of the power series 
expansion. 

In the past, the repeated total differentiation of a differential equation was con- 
sidered unfeasible, since, with increasing order, the derivatives become 
rather unwieldy expressions. Today, with the advent of fast electronic com- 
puters, such a procedure no longer seems necessarily unfeasible. It is well- 
known that in the last few years considerable progress has been made in auto- 
matic formula differentiation by computers. Apart from a straightforward 
differentiation of the differential equations, a great number of differential 
equations can be differentiated in a rather simple way after transforming 
them — by introducing auxiliary functions — into algebraic differential equations 
of the second degree. For special differential equations the procedure has 
been outlined in earlier papers by J. F. STEFFENSEN [20], E. RABE [17], 
and the author [10]. The procedure is based on the fact that the consecutive 
derivatives of a second-degree system of differential equations can be con- 
veniently obtained on a computer by recurrence formulas. 

12. The procedure is best illustrated by a simple example. Let us consider the 
differential equation 
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We introduce the auxiliary function 

e -x = u (2) 

and obtain from (1) and (2) a system of second-degree algebraic differential 
equations 


dx du 

dt U ’ dt 


-u 


Substituting the power series expansions 


(3) 


x = Z X v (t -^’ u= £ 


U v (t-to) V 


(4) 


v=o 


v=o 


n 


into (3) and comparing coefficients for the terms with (t - 1^,) results in the 
following recurrence formulas for the coefficients in (4): 


(n+l)X . = U 

v ' n+1 n 


n 

<” +1 > u „n = -Z U v U n^ 
v=0 


^ (n = 0, 1, 2, ...) 


(5) 


Since the first coefficient Xq is known from the initial value x(t^) for the step 
and the first coefficient U Q can be obtained from (2) , all following coefficients 
X v , U (v = 1, 2, 3, . . .) can easily be computed from the recurrence formulas 
(5) — a very convenient procedure for electronic computers. 


13. It is quite obvious that the power series expansion method allows — in an 
extremely simple way — an automatic step-size control. Assuming that we 
truncate the expansion (4) for x after the term X (t-1^) n , the leading term of the 
truncation error of x can easily be found by extending the computation to the next 
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coefficient X 


. If the truncation error turns out to be too large or too 

n+1 

small, the step size at At can immediately be adjusted in such a way that 

„ i 1 

lx (At) l remains within prescribed limits. For safety reasons it 
1 n+1' 

mi ght sometimes be advisable to consider more than just one term of the 
truncation error. In contrast to RUNGE-KUTTA or multistep methods, no 
repetition of any computation is necessary if the step size fails to meet the 
requirements for the magnitude of the truncation error. We know of no 
other method that offers such easy step-size control. 

14. Naturally, in our simple example No. 12, there is no real need to resort 
to the introduction of auxiliary functions, since a repeated differentiation 
of the differential equation (1) can be performed without difficulty. In Section 
. IV we shall present more involved examples that do not allow a convenient 
repeated straightforward differentiation without the introduction of auxiliary 
functions. 
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SECTION III. RUNGE-KUTTA TRANSFORMATION TYPE FORMULAS 


15. In two earlier papers [ll],[l2], the author presented RUNGE-KUTTA type 
formulas of high-order accuracy for the numerical integration of systems 
of first- and second-order differential equations. These formulas require 
a repeated total differentiation of the differential equations with respect to 
the independent variable. After m total derivatives are determined at the 
initial point t = ^, of the step under consideration, using , for instance, the 
method of Section n, a transformation of the dependent variables of the 
differential equations is performed in such a way that, in the case of second- 
order differential equations, the first m+2 total derivatives of these trans- 
formed dependent variables become zero for t = 

In the following we shall consider systems of second-order differential 
equations only, since these are the ones most frequently encountered in 
physics and mechanics. Moreover, our method is somewhat simpler in 
the case of second-order differential equations, because the number of 
RUNGE-KUTTA evaluations (including approximation of the truncation error) 
is reduced by 1 compared with the corresponding procedure for first-order 
differential equations. 

16. Let x be the original dependent variable — for the sake of brevity we shall 
write our formulas for one second-order differential equation 

x = f(t, x, x) (6) 

only, although they hold in the same way for systems . Let x T be the trans- 
formed dependent variable. Obviously the first m+2 total derivatives of x^ 
are zero for t = if we define 
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(7) 


m+2 

* T = x-^X v (t-t <5 ) V 
v=i 


with the power series coefficients X^ being defined as the v-th derivatives 
of x at t = t 0 , divided by the factorials v ! . 

From (7) it follows that 

m+2 

x T = x-^vx v (t-t 0 ) v_1 (8) 

V=1 


m+2 


x T = x-^v(v-i)X v (t-1 0 ) V ' 2 


v=2 

m+2 


•x T = f v(V-i)X v (t-t 0 ) V “ 2 = 


v=2 


Equation (9) represents the transformed differential equation 


x T = f T (t, x T , x T ) 


( 9 ) 


( 10 ) 


for which we derived, in papers [11], [12], RUNGE-KUTTA formulas of 
order m+4 as well as m+5. 

These formulas require no more than four RUNGE-KUTTA evaluations of the 
differential equations, including the determination of the leading term of the 
truncation error for x T . 

The small number of evaluations required is strictly a consequence of the 
fact that the first m+2 derivatives of x T are zero for t = to, since this be- 
havior of the derivatives drastically reduces the number of equations of 
condition for the RUNGE-KUTTA coefficients. 
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Since we have given a rather detailed derivation of our RUNGE-KUTTA for- 
mulas in papers [11] and [12], we restrict ourselves here to stating these 
formulas and to explaining, with the help of a flow chart, how to program 
our formulas on an electronic computer. 

17. First, we shall state the formulas derived in paper [11] for the case of second- 
order differential equations. The method described in this paper requires 
three RUNGE-KUTTA evaluations of (10) to obtain (m+4)-th-order formulas for 
x T and x T and one additional evaluation for an (m+5)-th-order formula for x r ^, 
m always standing for the number of differentiations performed on the original 
differential equation (6) at the initial t-value t = to of the current step. 

Using the traditional notation, the formulas of paper [ll] read: 


and 


k i = f T 0o + h, Xq, 0)h 
k 3 ^h, Xq + ftjkjh, 0 + pk x )h 

ks = f T (to + h, Xq, 0 + yk x + &k 2 )h 

^4 = frp(to + a 4 k , Xq + e 0 k x h, 0 + eki + £1^ + T]kg)h (h = step size) 


x T - Xq = Caok 2 h + 0(h B+s ) 

X rp - N} — (Cgokg "t Caok 3 + C^ol^h + 0(h B s ) > 
x - 0 — C 2 kg + C 3 k 3 + 0(h ) 


( 12 ) 


From (7) and (8) there follow for the initial valifes of the current step: 


x T (<o) = x(to)(=Xo), x T (to) = 0 (13) 

These values (13) have already been inserted into (11), 
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The first two equations in (12) yield two values, x T and x T , for t - t^+h 
that differ in accuracy by one h-power. Therefore, their difference can be 
considered an approximation of the leading term of the truncation error of 


The constants in formulas (11) and (12) are given by the following relatively 
simple expressions: 


m+2 _ m+1 

02 m+4 ’ 4 m+4 


2 /m+2\ m+1 = 3 m-2 /nrH \ m+1 

P ° ~ (m+4)® \m+4/ ’ 0 2 (m+2) (m+4) 2 Vm+4/ 



1 m-5, /m+l \ m+1 = 5 1 /m+lY 1 * 1 

2* (m+2) (m+4) \m+4y ’ 4 m+2 \m+2/ 




3 1 

— • 

4 m+4 





2 1 
3 (m+3)(m+5) 




It should be pointed out that the leading term of the truncation errors of x T 
and x T is not particularly small in our formulas. The situation is somewhat 
similar to that for the standard 4-th-order RUNGE-KUTTA formulas. In both 
methods , a part of the respective (m+5)-th-order or 5-th-order terms in the 
Taylor expansion for the solution is not covered at all. However, the extent 
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to which these leading terms of the truncation error are partially covered is 
very essential for the accuracy of such a method. A good coverage of these 
terms pays off by allowing a larger integration step size than in the case of 
a formula of the same order but with a larger truncation error. 

18. In paper [12] we presented formulas with smaller truncation errors than 
those in paper [11]. In fact, in the formulas in paper [i2], a parameter 
o is still available, and for sufficiently small values of a the absolute 
values of all members of the leading truncation error term, for x as well 
as for x , can be made as small as desired. However, they cannot be made 
zero, since some of the weight coefficients of our formulas would become 
infinite for a ^0. 

These more accurate formulas read as follows: 


k i = f T (to + K *o, 0)h 
ks = f T (to+ ^h, Xq+ fekih, 0 + 8k x )h 
k 3 = ^(to+^h, Xq+ Yo k ih + 6 0 kgh, 0 + yk x + 6kg) h 
K = f T (<o+ h, Xq, 0 + ek 1 + £k s + rik 3 )h 
and 

x T - Xq = (Caok 2 + Caok 3 )h + 0(h B+s ) 

x T - Xq = (<^+C 30 k 3 + c^h + 0(h ,+s ) i 

x T - 0 = C s k3+ C 3 k 3 +C 4 k4+0(h“ +5 ) 


(15) 


(16) 


Again, the initial values (13) have already been inserted into (15). 

The constants in (15) and (16) are no longer as simple as the expressions 
(14). Expressed as functions of the parameter a, these coefficients read: 
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m+2 

= 03 = 


1 -a 

1 + ^+ 4 ) 0 ® 


+ ? = <V“i§ 


v = _l._j_QL^ig 1 + (”? + 41 < ir2+mCT 1 + ( m+ ^ .' 

° 2 m+2 03 1 + (m+4)a 2 L 1 + (nrMjo 2 j 


Y = -°3 


m+i l~ q 2-(m+3)(m+4 
! m+2 [l + (m+4)o 2 ]® 


m+2\a 2 y 


/a s \” +1 1 + (m+4)a „ _ /a 3 \* +1 1 -a 3 - (m+2)(m+4)o 2 

\aj * a -( 1_CT ) [1 + ( m +4)o s ] 2 ’ \aj *m+2* [l + (m+4)cs 2 ] 2 


e = - 


m+2 1 - 2(m+4)a - (m+2)(m+4)d 2 


1.1.1 [3 + (m+2)o][3 - (2m + 13)a + 3(m+4)o a + (m+2)(m+4)(m+6)Q 3 ] 

£ “ a 3 m+1 *m+2"l - o' [3 - (m+2)(m+4)o 2 ][l - 2(m+4)a - (m+2)(m+4)o 2 ] 


1 cr 


[3 + (m+2)a][l + (m+4)a][l + (m+4)o 2 ] 


^ a 3 B+1 1-CT [3 - (m+2)(m+4)o 2 ][l -2(m+4)G- (m+2)(m+4)o 2 ] 

_ 1 m+5 1 2 - (m+4)Q - (m+2)(m+4)Q 8 

C 20 ~ ( m +2)(m+3)(m+4) ’ 1 - o' 3 - (m+2) (m+4)o® 

= _1 1 1 [1 + (m+4)q][l + (m+4)0 2 ] 

Cao a 3 “ +1 ’ (m+3)(m+4) *1 - a * 3 - (m+2)(m+4)o® 


C2 ~m+i 


C, = 


1 (m+5) s _ 1 < 2 - (m+4)0- (m+2)(m+4)o 2 

ag m+1 ’ (m+2) (m+3) (m+4) ' 1 - a ’ [3 - (m+2)(m+4)o 2 ][3 + (m+2)a] 

1 1 1 [1 + (m+4) 0 s ] 2 

a t3 m+1 ’ (m+3)(m+4)’ g(1-g)* 3 - (m+2)(m+4)o 2 

1 1 1 - (m+4)a - (m+4)(3m + lO)^ 3 - (m+2)(m+4) 3 a 3 


(m+3) (m+4) a 


[1 + (m+4)cr][3 + (m+2)o] 


/> _ 2 m+5 1 3 - (m+5) a - (m+2) (m+4) o 2 

Cao = a 3 ffi+1 ‘ (m+2) (m+3) (m+4)' 1 - a* [3 + (m+2) a] [3 - (m+2) (m+4) a 2 ] 


Cao 


2 1 1 # [1 + (m+4)0 2 ] 3 

a 3 " ,+1 * (m+3) (m+4) ' 1 - a ’ [1 + (m+4)a][3 - (m+2)(m+4)o 2 ] 

1 1 -2 (m+4) a - (m+2) (m+4) a 8 

" (m+3) (m+4) * [1 + (m+4)cr][3 + (m+2)a] 


In paper [12] these coefficients (17) are tabulated to 24 digits— for m = 3 
through m = 8. The parameter a is always chosen in such a way that the 
absolute values for the critical coefficients C 3 and C 4 are about 1/2 or less. 
The advantage of these more accurate formulas over the earlier formulas 
(11), (12), (14) will show up clearly in the examples in Section VI. 

19. It might be helpful to illustrate the application of our RUNGE-KUTTA pro- 
cedure by means of a flow chart, as is customary for computer programs. 
As Figure 1, we present a flow chart for our more accurate formulas (15), 
(16). The flow chart for our earlier formulas (11), (12) is almost identical 
with the flow chart in Figure 1. 
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FIGURE 1. FLOW CHART 








After having computed in 0 the power series coefficients X n (n=2, 3, . . . ,m+2) 
for the original differential equation (6) by the method indicated in Section n, 
then in 2 we introduce the proper arguments t, x^, x for the RUNGE-KUTTA 
evaluations of the transformed differential equation (9). However, since we 
have to deal with the original differential equation (6) instead of (9) we proceed 
in 3 from and x^ to x and x and evaluate in 4 , for these arguments, 

the right-hand side of (6) . In | 5 ) we compute the right-hand side of (9) , since 
our RUNGE-KUTTA procedure holds for the transformed differential equation 
(9) only. By multiplying by h we obtain in 6 the k t -values (15) for the 
transformed differential equation (9). After having computer all four k t -values 
(i=l, 2,3,4), we determine in 7 the approximate truncation error for x^. 

If necessary, the step size now has to be adjusted (by halving or doubling) in 
such a way that the truncation error remains within pre-set tolerances . After 
the step size has been checked and found to be satisfactory, the final values 
x and x for the end of the step are computed in j9 from the first and the 
third equation of (16). In 10 , at last, the final values for the original 

variables x and x are computed for the end of the step, and we are ready for 
the next step. 

20. The computational work in 3 and 10 of our flow chart is somewhat facili- 
tated by the fact that certain of the time increments in our RUNGE-KUTTA 
formulas are equal, namely tj-to = t 4 -^> = h. This means that the sums of 3 
and 5 have to be computed three times only, and no new computation of the 
sums in 10 is necessary. 

F urthermor e , according to 2 and 3 , X]_ = x*. This means that for the 

computation of f 4 in j 4 | the part of f that depends on t and x but not on x 
could be taken over from the computation of f x . In some cases — for instance 
equations (24) or (26) — this might practically reduce the number of evaluations 
of f by 1, as it does exactly in the case where f does not depend on x at all. 
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SECTION IV. THE EQUATIONS OF MOTION FOR THE 
RESTRICTED PROBLEM OF FOUR (THREE) BODIES 


The restricted problem of four bodies is based on the following assumptions: 

A. All four bodies (sun, earth, moon, space vehicle) are considered 
to be point masses that move in the same space-fixed plane. 

B. The center of mass of the earth- moon system moves with constant 
angular velocity (ju in a circle around the sun. 

C . The earth and the moon move with constant angular velocity in 
circles around the center of mass of the earth-moon system. 

D. The fourth body, the space vehicle, is of infinitesimal mass. Then, 
although attracted by the gravitational forces of the sun, the earth, and 
the moon, it can be considered to exert no gravitational forces on these 
bodies. 

In the space-fixed plane (Figure 2) we consider two different coordinate sys- 
tems: 

A. The space-fixed C-system with the sun S as origin (C = ? + i'll). 

B. The body -fixed z-system with the center of mass C of the earth-moon 
system as the origin and the direction from the earth to the moon as 
the x-direction (z = x + iy) . 


T1 



FIGURE 2. SPACE -FIXED PLANE 



As is customary in the restricted problem of three bodies, we choose the 
mass of the earth-moon system as the mass unit, the distance from the 
earth E to the moon M as the unit of distance, and the time unit in such a 
way as to make the angular velocity of the earth-moon system around its 
center of mass C equal to 1. This implies that the gravitational constant in 
Newton's gravity law becomes equal to 1. 


23. Let pi , 1-pi, pi, and m be the masses of the sun S, the earth E, the moon 

s 

M, and the space vehicle P, respectively, and let 0, Cg, and be their 
respective coordinates in the space-fixed C-system. To derive the equations 
of motion of the space vehicle we start from the Lagrangian function obtained 
for the space vehicle under the assumptions of the restricted problem of four 
bodies. Obviously, in the space-fixed ^-system the Lagrangian function 
L = T - U for the space vehicle reads as follows: 


L 



(18) 


It is customary in the restricted problem of three bodies to study the motion of 
the space vehicle in the body-fixed rotating z-system. We shall use the same 
z-system for the restricted problem of four bodies. The following relations 
obviously hold between the coordinates in the ^-system and the z-system: 


£ 

£ 

£ 


E 

M 


z e‘ (t+a) + Re“ 


i(t+ca T> iu,t 
-pi e ' + Re 


i(t+a) „ iwt 
(l-pi)e ' ; + Re 


> 


(19) 
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It follows from (19) that 


k-f E l ~ ~ l z -( 1 “H)l; kl = Iz+Re 1 ^ 1 ^ t °^| 

and by differentiation 

t _ .... . i(t+a) lout 

£ - (z+iz)e ' '+ Riuue 


or 

|f| = |(z+iz)+iuuRe^^ W ^ a ^| 


(21) 


Introducing (20) and (21) into (18), we obtain for the Lagrangian function L in 
the body-fixed rotating z-system 


L = 


m 

2 


| ( z+iz) + + “bL_ 

1 z+p| ]z-(l-u)| 


m|a c 


+ - 


|z + Re i[ <“ , - 1 > t - a] l 


(22) 


24. The equations of motion for the space vehicle in the case of the restricted 

problem of four bodies are then obtained from the Lagrangian equations of the 
second kind 







(23) 


by inserting expression (22) for L. 

The insertion results in the following equations of motion for the space vehicle 
in the case of the restricted problem of four bodies, as can easily be verified: 


« = o,Hv4. ( ,?r ,/ x+ l a ,, x-p. X + Rcoscp 

. y x Rcoscp-U -I ^ [(x + R CO s cp) s +(y-R sin cp) 2 ] 3 A 

y = -2x +y - ur*R sin cp- - U ~^ X + R cos^y^R sin cp f?/* 


} (24) 
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using the abbreviations 


\i = 1 - |i; cp = (l-aj)t + a 


(25) 


25. 


Omitting from (24) the terms contributed by the sun, our equations (24) 
reduce to the equations of motion for the space vehicle in the case of the 
restricted problem of three bodies: 


* = + x - ^ / [(x+J) 8 +%' s " ^(xVf+y^ 


y = -2x + y - ^ [(x+M f + /p/5 - ^ [(x Vf+y 8 ]*/ 5 


'i 


> 


(26) 


Equations (26) yield a first integral that can be obtained by multiplying the 
first equation in (26) by x and the second equation by f, then adding both 
equations and integrating with respect to time t. The result is 


|b^ + + y 2 )]- [(x+w'+y'iv 






7F* 


Const = J (27) 


the so-called Jacobi integral. 
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SECTION V. THE RUNGE-KUTTA TRANSFORMATION METHOD 

APPLIED TO THE 

RESTRICTED PROBLEM OF FOUR (THREE) BODIES 


26. 


Since we have explained our RUNGE-KUTTA transformation method in detail 
in Section III, we can now restrict ourselves to the problem of reducing the 
equations of motion (24) to an algebraic system of the second degree and de- 
riving the recurrence formulas for the coefficients of the power series ex- 
pansion for the solution of this second-degree system. 

In a completely obvious and straightforward manner, we introduce the follow- 
ing eight auxiliary functions into (24): 


cos cp = a, sin cp = b 
(x+nf+y 2 = p 2 , (x-p/f+y 3 = cf 3 , 




(x+Ra) s +(y-Rb) s = r 2 




> 


(28) 


Including the equations and differential equations for these auxiliary functions, 
we obtain the following second-degree system instead of (24): 


x = 2^+x+uu 2 Ra - u (x+|J.) - v (x-jj. 7 ) - w (x+Ra) 
y = -2x + y-upRb- uy-vy -w(y-Rb) 
a = -b(l-cu) 
b = a(l-cc) 

p 8 = (x-Hi) a +y® > 

cf = (x-n'f+y 3 

r 2 = (x+Ra) 2 + (y-Rb) s 


(29) 


pu + 3up = 0 


qv + 3vq =f 0 

rw + 3w r = 0 / 
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One might, at first glance, get discouraged when faced with ten equations. 
However, since system (29) is an algebraic system of the second degree, 
its numerical integration by power series expansions is completely trivial 
and is performed in exactly the same way as in the case of our simple example 
(1) in Section II. 


27. Let us denote the coefficients of the power series expansion of the functions 
x, y, a, b, p, q, r, u, v, w by the capital letters X v , Y y , B v , P v , Q v , 

R v’ U v’ V v’ W v’ res P ectivel y* 

By introducing the power series expansions 


x 


£X<t-y , 

V=0 



w 


Z w v“-W 

v=0 


and their derivatives, if required, into (29) and comparing coefficients of 
equal powers — for instance, coefficients of (t-to)" -1 or of (t-t 0 ) n — one obtains 
in a completely straightforward way, the following recurrence formulas for 
the ten functions occurring in (29): 


n-1 

(n+l)nX n+ i = 2 “VX n - 1 +a,2RA n - 1 -Z ( Uv +V v +W v)X n - 1 -v 

v=0 


n-1 


-^n-l + ^' V n-r R I W v A n -l-v 


v=0 


n-1 

< n+1)nY n + l = - 2nX n +Y n-r“‘ !RB n-l-Z (U v + V W v> Y n-l-v 

V=0 


n-1 


+ R 


VWB 1 

L. v n ~ 1 - 


■V 


v=0 


(30) 


nA = — (1— uu)B 
n n-1 
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For n = l, all sums in (30) with the lower limit v = i and the upper limit v = n-1 
have to be omitted. 

When we start the computation for a certain step, we know, from the initial 

values for this step, the coefficients X , X, and Y , Y . From (28) we then 

o 1 o 1 

obtain the coefficients A,B,P,Q,R,U,V,W. Setting n = 1 , we 

oooooooo 

find from (30) the coefficients X , Y , A , B , P , Q , R , U , V , W , and 

ZZ11X1111 1 

we continue by repeating the evaluation of the recurrence formulas (30) for 

n = 2, n = 3, etc. After having obtained in this way the coefficients X , X , 

z o 

. . . , x m+ 2 and Y 2 ’ Y 3 ’ • • • > Y m + 2 ’ we are read y f° r our RUNGE-KUTTA 
transformation procedure as described in Section III. 
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In the case of the restricted problem of three bodies we start from equations 
(26) instead of equations (24). Equations (28) reduce, in this case, to four 
equations only. In a quite obvious way, (29) then reduces to six second-degree 
equations, and (30) reduces to six recurrence formulas. 
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SECTION VI. SOME NUMERICAL RESULTS 

FOR THE 

RESTRICTED PROBLEM OF THREE (AND FOUR) BODIES 

28. In this section we present some of the numerical results we have obtained by 
applying the power series method (Section II) and our two RUNGE-KUTTA 
transformation methods (Section HI) to the computation of a periodic orbit of 
the restricted problem of three bodies, and to the computation of the corres- 
ponding orbit— with the same initial conditions— for the restricted problem of 
four bodies. 

29. Figure 3 shows, for the case of the restricted problem of three bodies, this 
periodic orbit in the rotating coordinate system. 


t»4 



FIGURE 3. PERIODIC ORBIT IN THE ROTATING COORDINATE SYSTEM 
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The periodic orbit was used in previous papers by the author and has the 
following initial values: 

Xq = 1.2, *o = 0, y 0 = 0, y 0 = -1.04935 75098 30320 (m = (31) 

The initial value y Q in (31) was obtained by an interpolation procedure that 
varied y 0 iteratively until the orbit, after half a period, finally crossed the 
x-axis perpendicularly. To preserve sufficient accuracy, the computation of 
y 0 for the periodic orbit was performed in 20-digit arithmetic. 

The computations presented in all tables in this paper were executed on an 
IBM 7094 digital computer (16 digits). 

Tables I and II list the results that we have obtained, in the case of the re- 
stricted problem of three bodies, by the methods described in Sections II and 
III for one orbit (about one month actual time), and for 12 consecutive orbits 
(about one year actual time), starting from the initial values (31). The pro- 
gram came to an automatic stop when the orbit intersected the x-axis again 
after one or 12 complete orbits, respectively. The last point of intersection 
with the x-axis was obtained by continuously halving the step size for the last 
step until we missed the x-axis by less than a pre-set tolerance (10~ 17 ). 

Table I refers to eighth-order formulas (m = 4) and Table II to twelfth-order 
formulas (m = 8). 


TABLE I. RESTRICTED PROBLEM OF THREE BODIES, 

RESULTS FOR EIGHTH-ORDER FORMULAS 

No. of No. of Computer Running 


Method* Orbits 

Final x 

Final y 

Steps 

Time (min) 

PSE 

1 

1.20000 00000 00040 

-1.049357509830345 

1594 

0.35 

RKT 1 

1 

1.20000 00000 00103 

-1.049357509830421 

1120 

0.24 

RKT2 

1 

1.20000 00000 00038 

-1.049357509830366 

840 

0.19 

PSE 

12 

1.20000 00000 00031 

-1.049357509830440 

19134 

4.20 

RKT 1 

12 

1.20000 00000 00430 

-1.049357509831333 

13459 

2.78 

RKT 2 

12 

1.20000 0000000010 

-1.049357509830525 

10080 

2.11 


(See footnote at end of Table II. ) 
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TABLE II. RESTRICTED PROBLEM OF THREE BODIES, 
RESULTS FOR TWELFTH-ORDER FORMULAS 


Method* 

No. of 
Orbits 

Final x 

Final y 

No. of 
Steps 

Computer Running 
Time (min) 

PSE 

1 

1.19999 99999 99981 

-1.04935 75098 30303 

493 

0.21 

RKT 1 

1 

1.20000 0000000001 

-1.04935 7509830321 

389 

0.15 

RKT 2 

1 

1.20000 00000 00013 

-1.04935 75098 30332 

290 

0.13 

PSE 

12 

1.20000 00000 00071 

-1.049357509830531 

5896 

2.49 

RKT 1 

12 

1.19999 99999 99991 

-1.04935 75098 30373 

4740 

1.83 

RKT 2 

12 

1.20000 00000 00097 

-1.04935 7509830627 

3353 

1.35 


' PSE 


Power Series Expansion Method [10] 


* < RKT 1 


Runge-Kutta Transformation Method [11] 


RKT 2 = Runge-Kutta Transformation Method [12] 


All methods listed in Tables I and n were run with automatic step-size con- 
trol for every step. The step size was accepted if for this step size At — but 
not for double the step size 2 • At — the absolute value of the truncation errors 
T x and T y in x and y were smaller than Ix^lO -16 , jy o |‘10~ 16 , respectively, with 
Xq and y 0 standing for the initial values for x and y for the current step. 
Comparing the computer running time in both tables, it becomes evident that 
twelfth-order formulas require considerably less computer time than eighth- 
order formulas — not to mention the prohibitively slow fourth-order standard 
RUNGE-KUTTA formulas. In fact, our twelfth-order formulas require only 
about 60 to 65 percent of the time for the eighth-order formulas. This time 
saving for the twelfth-order formulas is a consequence of the smaller number 
of steps required for a twelfth-order formula (only about 1/3 of the number 
required for an eighth-order formula in our example). 

Moreover, Tables I and II show a significant time saving for our RUNGE-KUTTA 
transformation method compared with the power series expansion method. Our 
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more recently developed RUNGE-KUTTA transformation method (RKT 2) gives 
results of about the same accuracy as the power series expansion method (PSE) 
in about half the computer time. 

As the x- and y-values of the tables show, all three methods are of about the 
same accuracy for equal tolerances; after 12 orbits we lost, in all methods, 
about 2 to 3 digits in x and y. This means that even after 12 orbits we miss 
the initial value Xq of (31) by only about 1/100 millimeter actual distance — 
certainly a negligible deviation. 

30. We now turn to the restricted problem of four bodies. For this problem we com- 
puted an orbit with the same initial conditions (31) that we used for our periodic 
orbit (Figure 3) in the case of the restricted problem of three bodies . The re- 
sults of our computations, if compared with our previous computations, will give 
an indication of how the attractive force of the sun affects our orbit of Figure 3. 
Naturally, the periodicity of our orbit is lost if the influence of the sun is taken 
into account. However, the shape of our orbit remains approximately pre- 
served for a surprisingly long time. Since the differential equations (24) for the 
restricted problem of four bodies are more involved than the corresponding 
differential equations (26) for the restricted problem of three bodies, the numeri- 
cal integration of equations (24) naturally took longer on the computer — about 
twice as long as for equations (26). 

In Tables III and IV we list our results for the restricted problem of four 
bodies. We use the starting values (31) and set or = 0 in (25). These assump- 
tions mean that sun, earth, moon, and space vehicle are initially all located 
on one straight line. We also made machine runs for other values for a. How- 
ever, because these changes in the configuration of the bodies do not produce 
any essential changes in our results, we can restrict ourselves here to pre- 
senting only the case wtiere cc = 0. 
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TABLE m. RESTRICTED PROBLEM OF FOUR BODIES, 
RESULTS FOR EIGHTH-ORDER FORMULAS 


Method* 

No. of 
Orbits 

Final x 

Final y 

No. of Computer Running 
Steps Time (min) 

PSE 

1 

1.19033 91358 54073 

-1.024871921155096 

1578 

0.67 

RKT 1 

1 

1.19033 91358 54174 

-1.024871921155201 

1117 

0.58 

RKT 2 

1 

1.19033 91358 54125 

-1.02487 1921155176 

835 

0.43 

PSE 

12 

1.17279 6212174518 

-1.0178210547 51673 

18144 

7.69 

RKT 1 

12 

1.17279 6212175284 

-1.0178210547 52938 

12954 

6.61 

RKT 2 

12 

1.17279 6212174671 

-1.0178210547 52064 

9723 

5.00 


* See footnote at end of Table n. 

TABLE IV. RESTRICTED PROBLEM OF FOUR BODIES, 

RESULTS FOR TWELFTH-ORDER FORMULAS 

No. of No. of Computer Running 

Method * Orbits Final x Final y Steps Time (min) 

PSE 1 1.19033 91358 54114 -1.024871921155165 485 0.38 

RKT 1 1 1.19033 91358 54016 -1.02487 1921155062 401 0.34 

RKT 2 1 1.19033 91358 54004 -1.024871921155048 287 0.27 

PSE 12 1.17279 6212174782 -1.01782 10547 52114 5598 4.35 

RKT 1 12 1.17279 6212174791 -1.0178210547 51654 5125 4.14 

RKT 2 12 1.172796212174737 -1.0178210547 51634 3280 2.73 

* See footnote at end of Table n. 


In our machine programs we again applied exactly the same automatic step-size 
control procedure, including the same tolerances, and the same procedure for 
the last (closing) step as in the case of the restricted problem of three bodies. 
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Since, in the case of the restricted problem of four bodies, the numerically 
correct values for our problem are not known, we can try to determine the 
accuracy of our methods only by comparing the results for the different methods. 
Lacking a better criterion for accuracy, it seems reasonable to assume that 
those digits in the final values for x and y which are in agreement for all 
methods are correct. But again, then, we do not lose more than 3 digits 
even after 12 orbits. This means, again, that the deviations in x among the 
different methods and the true solutions are, after 12 orbits, still of the order 
of l/l 00 millimeter. With respect to saving computer time — whether using 
twelfth-order formulas instead of eighth-order formulas or using our RUNGE- 
KUTTA transformation method instead of the power series method — we 
obtained about the same results as for the case of the restricted problem of 
three bodies. 

31. The orbit that we have considered is not perturbed very much by the influence of 
the sun, at least not for the first year for which we have run our computations. 
The deviations in x from Xq = 1.2 to x = 1.17279. . . after 12 orbits correspond 
to a deviation of about 10, 460 kilometers, which is less than one earth diameter. 
We also determined the x-deviations for the 1st, 2nd, 3rd, . . . , 11th orbit; 
they never exceed one earth diameter and seem to have an oscillatory behavior. 
However, periodic orbits of the restricted problem of three bodies which come 
closer to the moon than does our orbit turn out to be more sensitive to the 
influence of the sun and the moon. For such orbits, our model of the restric- 
ted problem of four bodies might no longer be sufficiently realistic. One 
might have to include the ellipticity of the moon orbit to obtain a satisfactory 
approximation of the real conditions. But, since this paper is mainly concerned 
with numerical integration procedures, we did not proceed further in the direction 
of a more realistic model. 
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The methods described in this paper are, however, applicable to considerably 
more involved problems than the restricted problem of three (or four) bodies. 
Actually, these methods have been applied in our Computation Laboratory to 
the problem of N oblate bodies as well as to the problem of the powered flight 
of a space vehicle — in both cases with rather favorable results. 


George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 
Huntsville, Alabama, June 1, 1966 
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APPENDICES 


32. In Appendices A through C we present summaries of some other numerical 
integration methods that we found interesting or promising and that are 
applicable to problems such as we have described in Section IV. 

Since these methods were developed and published by other authors, we 
restrict ourselves to presenting short descriptions without any derivation 
of the formulas. However, we shall give sufficient references to the origi- 
nal papers. 

In these appendices we also present, for comparison, some numerical re- 
sults obtained for these methods, applied to the problems of Section IV. 

We have tried to give an unbiased review of the methods in question, but it 
should be understood that we have based our opinion of these methods on their 
practical applicability to problems in celestial mechanics, such as we have 
described in Section IV. 
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APPENDIX A 


HIGH-ORDER EXPLICIT RUNGE-KUTTA FORMULAS 


33. In this appendix we consider the explicit RUNGE-KUTTA formulas of 
E. B. SHANKS [l 83, [19]. These formulas represent a remarkable ex- 
tension of the traditional 4-th-order RUNGE-KUTTA formulas to higher- 
order formulas. Naturally, the number of evaluations per step of the 
differential equations increases with the order of the formulas. However, 
the increase is not so sharp as to make high-order formulas uneconomical 
for an electronic computer. On the contrary, since higher-order formulas 
permit — without loss of accuracy — a larger step size than 4-th-order for- 
mulas, the differential equations can be integrated in considerably fewer 
steps . This more than compensates for the increased computational 
effort per step in such high-order formulas . 

34. The most interesting and highest-order formula of E. B. SHANKS is an 8- 
th-order formula based oh 12 evaluations of the differential equations . We 
restrict ourselves to quoting this 8-th- order formula and applying it to the 
problems of Section IV. 

E. B. SHANKS’ 8-th-order formula reads, when written in the traditional 
notation: 
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K = 
K = 
K = 

k 4 = 

^5 = 

K 

k? 

K 

K = 

kjo = 
= 

^12 = 


f (t 0 ,Xo)h 

f ( t ° + f h ’ x ° + ^ k l) h 

t( t ° + |h.S» + 5ik 1 + ik=) h 

f(‘o + jh.% + ^ki + ^k 3 )h 


\ 


("° + W h - *o + H ^ + H ^ ^l h 


125 


X 0 + ^ k 1 + 5IJ k 4 + ^|l%) h 


243 "* 972 ; 

125 


= £(‘o + |h, 

= £(<b + |h.Xo-^ki + y k 4 + ^-ks-!l%) h 

■ 4 + f h '^-i! ki -il k ^ + ii^ + M k ’) h 
(^ + i h ' 


, 1175, 32, 3125, , , 121, 

x o + ^r k i-T k 4-T^r k e + 26k 6 + 7^ 


f (t,+| h >!!o +||| k 1 -|ik,- 
(t„ + |h, 


162 

1375 


162 

59 


i 4. 51, 

324 k® 9 ^ 162 


-H^) h 

k 7 + | k 8 + 1 %) 1 


, 1303 , 71, 1375 , 37 , A 103, . 1 . V 

^ 1620 kl 27 k4 324 ^ + 6 ^ 162^ + 10 


L , , 955, , 2560, , 8125, 612, , 7 , 

h>X ° " 492^ + ^" k 4 + ^^ k s-— k 6 + 77 k ? 


369 


738 


41 


162 

l-l 

82' 


27 

164 


ks-if‘%-f k B + f k u) h 


> (A-l) 


x = Xq+ — (41^+ 2161% + 2721^+ 27kg +27159+361^0 +180^+411^2) + 0(h 9 ) 


/ 


35. We have programmed SHANKS’ formulas (A-l) for the restricted problem of 
three — as well as four — bodies. Since no better step-size control procedure 
seems to exist for SHANKS' formulas, we have applied RICHARDSON’S de- 
ferred approach to the limit. This, naturally, is a considerable additional 
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computational effort merely for the benefit of step-size control. 

We have run the orbits of Section VI with SHANKS' formulas (A-l) on the 
same electronic computer and with the same tolerance 10 -16 that we used for 
the various methods for which we listed numerical results in Section VI. 

The following Tables, A-I and A-II, present the results of the runs we made 
with SHANKS' formulas for one orbit. 

TABLE A-I. RESTRICTED PROBLEM OF THREE BODIES 

No. of No. of Computer Running 

Method Orbits Final x Final y Steps Tim e (min) 

RKS* 1 1.20000 00000 00002 -1.04935 75098 30310 814 0.46 


TABLE A-H. RESTRICTED PROBLEM OF FOUR BODIES 

No. of No. Computer Running 

Method Orbits Final x Final y Steps Time (min) 

RKS* 1 1.19033 91358 54033 -1.024871921155064 817 1.63 

* RKS = RUNGE -KUTTA -SHANKS (8-th-order formula) 

The accuracy of SHANKS' formulas (A-l) is quite impressive. Comparison of 
the final values in Table A-I with the initial values (31) shows that we lose 
only 1 digit in x and 2 digits in y. This is somewhat less than we lose in 
Table I for the power series expansion method and our RUNGE-KUTTA trans- 
formation methods when set up as 8-th-order methods . The running time for 
SHANKS' formulas, however, is considerably longer than for the methods of 
Tables I and III. In the restricted problem of three bodies SHANKS' method 
takes about twice as long as our RKT 2 method, and in the restricted problem 
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of four bodies almost four times as long. The numerous evaluations — required 
by SHANKS' method — of the trigonometric functions sin cp and cos cp, occurring 
in equations (24), account for the relatively long running time for this method 
in the case of the restricted problem of four bodies. The presence of trans- 
cendental functions in the differential equations will always slow SHANKS' 
method, since in his 8— th-order formula these functions have to be evaluated 
23 times per step (including step-size control procedure) versus 4 evaluations 
per step for our RUNGE-KUTTA transformation formulas. 

Naturally, in our methods in Sections H and HI we must also pay for the com- 
putation of the derivatives that are required in these methods. However, the 
computation (especially of the lower-order derivatives) is rather easy and 
fast by the use of our recurrence formulas. 

SHANKS' formulas might gain considerably if a less expensive but still reliable 
step- size control procedure were available for them. 

However, since SHANKS' formulas are of 8-th order at best, one cannot ex- 
pect them to compete with our higher-order formulas, as a comparison of 
Tables A-I and A-II with the first part of Tables II and IV clearly indicates. 
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APPENDIX B 


HIGH-ORDER IMPLICIT RUNGE-KUTTA FORMULAS 


36. In explicit RUNGE-KUTTA formulas — for instance, Formula (A-l) — the com- 
putation of the increment k v requires a knowledge of only the preceding in- 
crements kj, kg, ... k v _ x . Therefore, in explicit RUNGE-KUTTA formulas, 
all increments k v can be computed one after the other in one procedure. 

In implicit RUNGE-KUTTA formulas the increment k v depends not only on the 
preceding increments 1%, kg, . . . k v-1 but also on k^ itself and on the succeed- 
ing increments k^, k^, .... Therefore, in implicit RUNGE-KUTTA for- 
mulas, the increments k^, have to be determined by an iterative procedure. 
Naturally, such an iterative computation is more involved than the straight- 
forward procedure for explicit RUNGE-KUTTA formulas. 

However, there are some points in favor of implicit RUNGE-KUTTA formulas. 
For instance, implicit RUNGE-KUTTA formulas are available for any order, 
whereas no explicit RUNGE -KUTTA formulas exceeding the 8-th order are 
known so far. 

37. Implicit RUNGE-KUTTA formulas have been studied by J. KUNTZMANN and 
his collaborators. We mention as a reference the textbook of F. CESCHINO 
and J. KUNTZMANN [9 ]. More recently, J. C. BUTCHER has published 
two noteworthy papers on implicit RUNGE-KUTTA methods. In these 
papers he derived various implicit formulas based on the quadrature formula 
of GAUSS-LEGENDRE [4] or on quadrature formulas of RADAU [5]. The 
latter formulas have the advantage of requiring fewer iterative kv stages than 
the formulas based on the GAUSS-LEGENDRE quadrature formula. In a 
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separate paper [7], J. C. BUTCHER presents 20-digit tables for the coeffi- 
cients of his implicit RUNGE-KUTTA formulas (up to the 20-th order). 


38. Let us illustrate the procedure for implicit RUNGE-KUTTA formulas by quot 
ing one of BUTCHER'S 8-th-order formulas based on RADAU's quadrature 
formulas: 


K = f(*o) h 

kg = f(XQ + PsiL + fegkg + (3g3k 3 + 034.1%) h 

k 3 = f(Xo+ 031 k x + 033^2+ 033 k 3 + 034^) h 

1^4 = f (Xq + 041 kj+ 04 gkg + 043 kg + 0441%) h 

ke = f(Xo+ 0 a k!+ $ 5 gkg+ 0531%+ 054l%)h 

x = 5%+ Ciki + C s kg+ C 3 kg+ C 4 1%+ C 5 k 5 + 0(h 9 ) 




> (B-l) 


/ 


It is assumed that the independent variable t does not appear explicitly on the 
right-hand side of the differential equation. This is no restriction, however, 
since by the introduction of an additional dependent variable (which is identical 
with t) the independent variable t can always be eliminated on the right-hand 
side of the differential equation. 

From (B-l) it follows that these implicit RUNGE-KUTTA formulas, which 
correspond to SHANKS' formulas (A-l), have only three iterative stages (kg, 

1%, 1%) and only five stages altogether. 

Formulas of 12-th-order of the type (B-l) would contain five iterative stages 
and seven stages altogether. 

In (B-l) one first computes k x and starts the iteration for k 2 , k 3 , 1% by introduc- 
ing into the right-hand side of the 2nd, 3rd, and 4th equation the value k x as the 
first approximation of kg, k 3 , 1%. After the iteration procedure has converged 
to the final values kg, k 3 , 1%, the last two equations of (B-l) yield ks and the 
x-value for the end of the step. 
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J. C. BUTCHER has pointed out that a step-size control procedure is possible 
for formulas of the type (B-l) by combining these formulas with implicit 
RUNGE-KUTTA formulas of the same order but based on the GAUSS- 
LEGENDRE quadrature formula. This means, however, that the computa- 
tional effort has to be doubled to obtain a reliable step-size control. In this 
respect, one is faced with the same situation as in the case of explicit 
RUNGE-KUTTA formulas. 

39. Using BUTCHER'S formulas of the type (B-l), and for step-size control the 
corresponding formulas of the GAUSS-LEGENDRE type, we again computed 
the orbit of Section VI for the restricted problem of three bodies. 

Table B-I presents the results of the runs we made with BUTCHER'S 8-th- 
and 12-th-order formulas. 

TABLE B-I. RESTRICTED PROBLEM OF THREE BODIES 


Method 

No. of 
Orbits 

Final x 

No. of Computer Running 
Final y Steps Time (min) 

RKB (8)* 

1 

1.20000 00000 00010 

-1.04935 75098 30318 870 

1.56 

RKB (12)* 

1 

1.20000 00000 00013 

-1.049357509830328 216 

0.88 


RKB (8) = RUNGE-KUTTA-BUTCHER (8-th-order formula) 

RKB (12) = RUNGE-KUTTA-BUTCHER (12-th-order formula) 


The runs were made on the same computer and with the same tolerances for 
the truncation error as the runs reported in Section VI and in Appendix A. The 
values in Table B-I can be compared with the first half of Tables I, H, and 
A-I. Quite obviously, in our example, the implicit RUNGE-KUTTA formulas 
are much slower than the corresponding explicit RUN GE -KUT TA formulas or 
the RUNGE-KUTTA transformation formulas. 
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However, since implicit RUNGE -KUTTA formulas are available for any order 
of accuracy, they might be rather attractive for some problems which require 
high-order accuracy. 
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APPENDIX C 


MODIFIED MULTISTEP METHODS 


40. Although this paper deals with one-step methods, for comparison we shall 
consider in this appendix the new modified multistep formulas suggested 
by J. C. BUTCHER [8 ]. We have already mentioned in Section I that 
these new formulas are much more convenient than the traditional multi- 
step formulas since they are based on fewer backward steps than were the 
earlier formulas. 

For instance, BUTCHER'S 7-th-order formula is based on only three equi- 
distant backward points t„, t n _ 1; ta-g, whereas the traditional 7-th-order 
implicit multistep formula of ADAMS requires six equidistant backward 
points t n , t,^, t n _ s , t n _ 3 , t n _ 4 , t^g. 

As already pointed out in Section I, the stability restrictions of the traditional 
multi step methods are overcome in the case of the modified multi step methods 
by the introduction of one additional non- step point. This requires one addi- 
tional formula for the non-step point. But this additional computational effort 
seems more than compensated for by the convenience of the new formulas. 
Since they are based on fewer backward points, the starting procedure and 
the change of interval size is much more easily performed for these formulas 
than for the traditional multistep formulas. 

41. BUTCHER' s multistep formulas are of the predictor-corrector type consisting 
of two predictor formulas for the non- step point and for the next step point, and 
one corrector formula for the next step point. The formulas require three 
evaluations of the differential equations per step. In his paper [8 ], BUTCHER 
presents such multistep formulas of the 5-th, 7-th, 9-th, 11-th, and 13-th 
order. 
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For instance, let us consider his 7-th-order formulas. For a first-order 
differential equation x = u(t, x) his formulas read: 
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The subscript n + l/2 denotes the non-step auxiliary point that is required 
for stability reasons. 

These formulas (C-l) were also used in an iterative way when starting the 
method and when halving the integration step size. 

In the latter case, we proceeded from the 5-th-order formula 
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as a first approximation for the half step point. 

In the case of a second-order differential equation x = f(t, x, u), a convenient 
approximation for the leading term of the truncation error for x n _ L m (C-l) 

reads: 
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The formula (C-3) is obtained by constructing an 8-th-order formula for 
x n+1 which coincides with the third formula in (C-l), as far as the x-terms 
are concerned 


135 . 31 , 

x n+x = x n + — (x n - x n _!) + — (x a - x n _a) 


617 


+ 


•as- 

■(- 


405 27 1 \ 

‘ - - ~ -< ao >i — 1 -» no/t ^-n - §] 


+1 1234 “ 1234 n_1 1234 


(C-4) 


+ hS/ 51 , . 351 , . 81 


-f + ^-f + 

, 1-n 4.1 ~ . in 1 


1234 n+1 1234 n 1234 

and subtracting (C-4) from the third equation of (C-l). 
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42. We have applied BUTCHER’S 7-th-order formulas to the restricted problem 
of three bodies. The point of intersection with the x-axis of our periodic 
orbit (Figure 3) after a complete period was found from the coordinates, 
velocities, and accelerations for the last three points of the orbit by 
Hermitian interpolation. 

While all results in Section VI and in Appendices A and B were obtained using 
an error tolerance of 10~ 16 , we had to relax the tolerance somewhat in the 
case of BUTCHER'S multistep method. In view of the magnitude of some of 
the coefficients of the predictor formulas, the necessity for such a relaxa- 
tion is not surprising. We ran BUTCHER'S formulas for an error tolerance 
of 5‘10~ 16 and of 1(T 15 . 
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TABLE C-I. RESTRICTED PROBLEM OF THREE BODIES 


Results for 7-th-Order Formulas 
Tolerance: 5 • 10" 16 


Method 

Final x 

Final £ 

No. of 
Steps 

Computer Running 
Time (min) 

BUTCHER 

1.19999 99999 99904 

-1.04935 75098 31055 

3370 

0.63 

PSE 

1.20000 00000 01013 

-1.04935 75098 31284 

2118 

0.38 

RKT 2 

1.20000 00000 00124 

-1.04935 75098 30410 

1097 

0.20 


TABLE C-H. RESTRICTED PROBLEM OF THREE BODIES 

Results for 7-th-Order Formulas 
Tolerance: 10“ 15 

No. of Computer Running 


Method 

Final x 

Final y 

Steps 

Time (min) 

BUTCHER 

1.20000 00000 01204 

-1.04935 75098 31504 

2799 

0.53 

PSE 

1.200000000002329 

-1.049357509832640 

1913 

0.35 

RKT 2 

1.20000 00000 00239 

-1.04935 75098 30579 

989 

0.18 


Tables C-I and C-II show the results of the runs we made for BUTCHER'S 
modified multistep formulas and for certain other methods described in 
this report, always comparing methods of the same order and runs with 
the same error tolerance. 

Since a multistep method is based on a number of backward steps, it neces- 
sarily has a larger truncation error than a one-step method. A multistep 
method therefore requires a larger number of steps for our orbit. 

In Tables C-I and C-II the computer running times are rather closely propor- 
tional to the number of steps, independently of the individual method. 
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For BUTCHER'S multistep method this means that, although this method re- 
quires only three evaluations of the differential equations per step, it is not 
faster per step than the one-step methods we have considered. Obviously, 
much of the speed of the multistep method is lost by the frequent step size 
changes required in our problem. Any halving of the step size requires an 
expensive (forward and backward) iteration procedure, to maintain sufficient 
accuracy. Therefore, multistep methods are not well-suited to problems 
that require frequent changes in the integration step size. Under the circum- 
stances it is still somewhat surprising how well BUTCHER'S multistep 
method is doing in our example compared with the one-step methods listed. 
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